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I.  INTRODUCTION 

Eddy-current  testing  was  first  systematically  studied  in  Germany 
during  World  War  II,  but  did  not  receive  wide  recognition  until  Forster 
and  his  colleagues  published  the  results  of  their  extensive  theoretical, 
experimental  and  industrial  investigations  between  1952  and  1954  [1]. 
These  papers  did  not  include  a  quantitative  theory  for  flaw  detection, 
however,  and  it  was  not  until  1964  that  Burrows  constructed  such  a 
theory  in  his  dissertation  [2].  His  idea,  which  was  based  on  electro¬ 
magnetic  theory,  is  to  replace  the  flaw  by  an  equivalent  dipole,  and  is 
reasonable  if  the  flaw  is  small  compared  with  the  skin  depth  (sometimes 
called  the  diffusion  length)  of  the  eddy-currents. 

The  work  of  Dodd,  et  al.  [3-5]  brought  eddy-current  nondestructive 
evaluation  (NDE)  "of  age"  by  showing  how  one  could  get  useful  analytical 
results  based  on  a  rigorous  application  of  electromagnetic  theory. 

Their  theory  of  flaw  detection  contains  Burrows’  and  is  subject  to  the 
same  limitations  of  dipole  representation.  Equally  significant  with 
their  use  of  rigorous  electromagnetic  theory,  in  our  opinion,  is  their 
development  of  computer  programs  to  compute  the  integral  representations 
of  the  electromagnetic  fields  [6-7].  Now  the  modem  era  of  eddy-current 
inspection  is  upon  us,  based  on  the  union  of  numerical  methods  in 
mathematics  and  rigorous  electromagnetic  theory. 

As  further  evidence  of  this,  we  cite  the  extensive  research  sup¬ 
ported  by  the  Electric  Power  Research  Institute  (EPRI)  [8],  in  par¬ 
ticular  the  work  of  W.  Lord,  of  Colorado  State  University,  B.  A.  Auld, 
of  Stanford  University,  and  A.  N.  Mucciardi,  of  Adaptronics,  Inc. 

Lord  is  developing  a  finite  element  model  for  eddy-current  nondestruc¬ 
tive  testing  phenomena,  whereas  Auld  is  using  the  reciprocity  theorem 
of  electromagnetics  to  quantitatively  model  flaw  responses  in  eddy- 
current  testing.  Finally,  Mucciardi  is  developing  a  system  for  flaw 
detection  and  classification  by  using  an  adaptive  learning  network  for 
eddy-current  signal  analysis. 

In  this  report  we  describe  an  approach  to  the  reconstruction  of 
flaws,  not  merely  their  detection.  This  will  give  us  the  ability  to 
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obtain  much  more  information  about  the  nature  of  the  flaw,  unimpeded  by 
the  size  restriction  of  the  dipolar  approximation  that  was  mentioned 
above.  By  "flaw"  we  mean  virtually  any  departure  of  the  medium  from  a 
standard  condition,  which  is  known  a  priori,  such  as  may  be  produced  not 
only  by  a  crack  but  also  by  conductivity  in  homogeneities  produced  by 
stresses,  magnetite  build-up,  etc.  Our  approach  is  very  much  in  the 
spirit  of  contemporary  work  in  invevse  methods  in  electromagnetics  [Si¬ 
ll]  and  electromagnetic-geophysical  prospecting  [12-19]. 

At  this  point  we  introduce  some  systems-related  ideas  that  should 
make  clearer  the  way  our  concept  of  inversion  is  to  be  used  for  non¬ 
destructive  evaluation.  Refer  to  Figure  1,  which  shows  a  "system," 
together  with  its  input  and  output.  In  part  (a)  of  the  figure  the  input 
is  known  and  so  is  the  system,  and  the  output  is  to  be  determined.  This 
is  the  "forward"  or  "direct"  problem.  For  example,  the  input  could  be  a 
current  or  voltage  source,  and  the  system,  a  coil  coupled  to  the  work- 
piece.  The  output,  the  magnetic  vector  potential  or  induced  eddy- 
current  within  the  workpiece,  can  be  directly  computed  in  a  straight¬ 
forward  manner  by  electromagnetic  theory. 

In  part  (b)  the  system  and  output  are  known,  ani  the  input  is  to 
be  determined.  This  is  a  problem  of  communication  theory,  or  signal 
detection.  One  assumes  a  catalog  of  possible  input  signals  to  be  avail¬ 
able,  whose  structure  and  characteristics  are  known  a  priori;  from  the 
known  output  one  estimates  the  input  signal  on  the  basis  of  the  maximum 
a  posteriori  probability  of  its  occurrence.  This  is  the  basis  of  Adap- 
tronics’  flaw  detection  system.  Their  inpi’t  data  base  (the  catalog  of 
possible  input  signals)  consists  of  506  defect  waveforms  under  tube 
supports  and  261  isolated  defect  waveforms,  all  of  which  are  at  400  kHz. 

This  is  an  example  of  another  "forward"  method,  and  appears  to  be 
sufficient  for  many  applications.  It  is,  however,  limited  by  both  the 
large  volume  of  signal  waveforms  that  must  be  catalogued  for  a  suitable 
gereralized  interpretation  and  by  the  subjective  comparisons  made  by  the 
interpreter.  The  method  also  gives  little  indication  of  the  sensitivity 
of  the  solution  to  possible  errors  in  the  data  and  the  degree  of  non¬ 
uniqueness  associated  with  the  chosen  mode. 


2 


In  Figure  1(c)  both  the  input  and  output  are  known  and  the  system 
is  unknown.  The  input  could  be  a  known  "probing"  signal,  and  the  out¬ 
put,  the  measured  response  to  the  probe.  The  object  is  to  determine  the 
nature  of  the  system.  This  is  an  example  of  "system  identification"  or 
"parameter  estimation,"  where  "parameter"  refers  to  certain  parameters 
of  the  unkown  system.  In  the  sense  that  problems  (a)  and  (b)  are 
"direct,"  problem  (c)  is  the  "indirect"  or  "inverse"  problem,  and  is  the 
problem  discussed  in  this  report.  What  one  seeks  in  this  problem  is  a 
model-system  that,  when  operating  on  the  input,  produces  a  "modex- 
output"  that  is,  in  some  sense,  an  optimum  estimation  of  the  known  output. 

There  are  two  common  ways  of  obtaining  a  sufficient  amount  of  inde¬ 
pendent  data  to  estimate  parameters  in  eddy-current  testing:  (a)  through 
the  use  of  multiple  coils,  and  (b)  through  the  use  of  multiple  fre¬ 
quencies,  including  pulses  or  transient  signals  [20],  Of  course  a  com¬ 
bination  of  the  two  may  be  used.  In  this  report  we  consider  only  the 
multiple  coil  method. 

In  Figure  2  we  show  a  system  of  coaxial  coils  within  a  tube  that  is 
to  be  inspected.  Within  the  wall  of  the  tube  is  located  an  anomalous 
region  (the  "flaw")  that  we  wish  to  reconstruct.  A  mathematical  mesh 
is  defined  that  surrounds  the  anomaly,  as  shown.  The  properties  of  the 
mesh,  such  as  its  location,  size,  and  fineness,  are  known  to  us.  Wnat 
we  don't  know  are  the  values  of  the  electrical  conductivity  to  assign  to 
each  rectangle  of  the  mesh.  The  "system,"  then,  consists  of  the  coils, 
the  tube,  and  the  mesh  that  encloses  the  anomalous  region,  .‘s  unknown 
parameter^  that  are  to  be  estimated  in  order  that  the  system  be  "iden¬ 
tified,"  in  the  sense  of  Figure  1(c)  and  its  discussion,  are  the  conduc¬ 
tivities  that  are  to  be  assigned  to  each  rectangle  of  the  mesh.  The 
known  input  is  the  current  to  the  exciting  coil,  and  the  known  outputs 
are  the  voltages  induced  into  the  sensing  coils.  Clearly,  if  we  can 
determine  the  conductivity  map  that  is  defined  on  the  mesh,  we  will  have 
reconstructed  the  anomalous  region. 

The  method  of  solving  this  problem  is  based  on  minimizing  the 
square  of  the  error  between  the  actual  measured  data  and  that  produced 
by  the  model-system,  the  model-output  (this  error  is  often  called  the 
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residual) .  The  parameters  that  are  varied  to  produce  the  optimum  model, 
in  the  least-squares  sense,  are,  of  course,  the  conductivities  that  are 
assigned  to  each  cell  in  the  mesh  of  Figure  2. 

Thus,  mathematically,  we  wish  to  determine  a  set  of  unknown  para¬ 


meters  0^,  j=l,  .  .  .  ,  M,  where  M  is  the  number  of  cells  in  the  mesh, 
from  a  set  of  data,  i=l,  .  .  .  ,  N,  where  e^  are  the  voltages  induced 


into  the N sensing  coils.  The  e^  are  functionally  related  to  the  a in 


a  known  way;  that  is 
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Hence,  given  the  wa  can  calculate  the  e^  by  treating  this  as  a 
"forward"  problem,  in  the  sense  of  Figure  1(a).  The  equations  (1)  that 


define  the  forward  problem  are  determined  by  using  electromagnetic  theory. 


But  it  is  the  voltages,  e. ,  that  are  the  given  data,  so  we  muse 


invert  the  system,  (1) ,  to  determine  the  .  We  do  this  by  minimizing 
the  error  function 


F(a1 


N 
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a  )  =  [  Z  (e  - 
M  i=l  1 
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Iterative  methods  are  commonly  used  to  carry  out  the  minimization 
of  (2).  The  iterative  method  successively  improves  a  current  model, 


i.e.,  a  current  estimate  of  the  a^,  until  the  error  measure,  (2),  is 


small  and  the  parameters  are  stable  with  respect  to  reasonable  changes 
in  the  model. 

The  success  of  this  method  of  inversion  depends  largely  on  the 
availability  of  suitable  numerical  algorithms  for  carrying  out  the 


least-squares  solution  of  (1).  Any  algorithm  chosen  m-st  contend  with 
the  fact  that  the  problem  as  posed  in  (1)  and  (2)  is  generally  quite 
ill-conditioned,  which  means  that  small  variations  in  input  data  can 
produce  quite  large  variations  in  the  solution.  The  commercially  avail¬ 
able  FORTRAN  packages,  UNPACK  (21]  and  MINPACK  [22],  contain  well- 
written  codes  for  least-squares  alogrithms,  and  these  codes  served  as 
the  basis  of  the  numerical  experiments  to  be  described  in  this  report. 
UNPACK  consists  of  linear  equation-solving  algorithms,  and  MINPACK 
contains  nonlinear  least-squares  algorithms.  A  description  of  these 
algorithms  will  be  given  3n  a  later  section  of  this  report. 

These  experiments  indicate  that  the  inversion  method  works  quite 
well  on  simulated  flaws,  even  when  the  data  is  corrupted  by  as  much  as 
20%;  this  is  quite  important  in  applications.  Another  nice  feature  is 
that  once  the  nonlinear  inversion  algorithm  has  converged,  it  is  possible, 
using  the  techniques  of  linear  inverse  theory,  to  assess  the  errors  and 
resolution  in  the  estimate  of  the  final  model.  The  objective  is  to 
determine  which  features  of  the  model  are  well-resolved  and  important 
to  the  interpretation  of  the  data  and  which  features  are  irrelevant, 
in  the  sense  that  the  data  neither  support  nor  reject  their  inclusion  in 
the  model.  This  is  also  quite  useful  in  eddy-current  NDE. 


II.  DERIVATION  OF  THE  MODEL  SYSTEM 


(a)  Integral  Equations 

The  detection  of  flaws,  or  anomalies,  by  means  of  eddy-currents  de¬ 
pends  upon  the  fact  that  flaws  are  not  electrically  conducting  and  that 
the  eddy-current  flow  is  interrupted  at  the  boundary  of  the  flaw.  The 
flaw,  therefore,  can  be  considered  to  be  an  inhomogeneity,  which  consists 
of  a  conductivity,  o^,  that  is  imbedded  in  a  region  whose  conductivity, 
Oq,  is  known  a  priori.  The  dielectric  constant  and  magnetic  permeability 
of  each  region  are  those  of  free  space,  and  y^.  Hence,  Maxwell's 
equations  for  the  two  regions  are: 
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V  X  Eg  =  -jcoPgHg  (Known  Region)  (3)  (a) 

V  x  H0  =  (aQ  +  j“£0)E0 

V  x  Ef  =  (Flawed  Region)  (3)  (b) 

V  x  af  =  (Cf  +  jwe0)Ef 

Upon  subtracting  (3)(b)  from  (3) (a),  we  get 

7  *  «0  -  v  -  i“V5o  -  V 

V  X  (H0  -  a£)  ■  a0( E0  -  Ef)  +  jue0(E0  -  Ef)  +  (o0  -  Of)Ef  .  (4) 

/e  have  added  nd  subtracted  O^E^  to  get  the  final  form. 

Thus,  the  perturbation  of  the  electromagnetic  field,  E^  -  Ef, 

5^  -  Hf,  satisfies  the  s»me  equation  as  the  original  electromagnetic 
field  within  the  known  region,  except  for  the  presence  of  the  terra 
(Og  -  af)E£.  This  term,  which  is  equivalent  to  a  current  source,  J^, 
iepresents  the  presence  of  the  anomalous  region,  or  flaw.  It  is  impor¬ 
tant  to  note  that  J  vanishes  off  of  _ne  flaw,  because  there  a,  = 

a  f_  0_ 

In  the  usual  way  we  can  derive  a  vector  wave  equation  for  E^  -  E^ 
from  (4) : 


V  x  V  x  (EQ  -  Ef)  =  -jajy07  x  (HQ  -  Hf) 


(a)  pq£0  -  jo)pnan)(Ert  -  Ef)  -  ju)Pn(or,  -  . 


0  O' v  0 


0V  0  f'  f 


(5) 


By  treating  the  last  term  in  (5)  as  a  source,  we  can  immediately 
write  down  a  formal  solution  fcr  the  perturbed  field  Eq  -  E^: 


ce0  -  e£><;>  -  jWo[| 


G(ijP)  .  If(?’)(o0  -  of)(i’)dP  ,  (6) 


Flaw 
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where  G(rjr’)  is  the  dyadic  Green's  function  for  the  known  region  [23], 
and  the  volume  of  integration  is  the  flawed  region,  for  which  0^  -  0^ 

#  0.  This  equation  is  the  basis  for  our  inversion  technique.  Before 
going  further  with  it,  we  make  the  following  observations  which  will 
allow  us  to  reduce  the  problem  to  a  scalar  system. 

In  eddy-current  NDE  work,  we  don't  measure  perturbed  fields  di¬ 
rectly;  rather,  uie  perturbed  EMF  that  is  induced  in  a  probe  coil  is 
measured.  Such  an  EMF  is  given  by  the  line  integral  of  (6)  along  the 
probe  coil  windings.  In  the  system  that  we  are  investigating,  as 
shown  in  Figure  2,  the  EMF  is  the  line  integral  of  the  azimuthal  elec¬ 
tric  field,  E^  (in  the  usual  cylindrical  coordinate  system,  wherein  the 
z-axis  coincides  with  the  axis  of  the  tube).  Hence,  all  we  need  to 
consider  is  that  single  component  of  the  E  vector.  In  addition,  we  can 
ignore  any  ^variations  of  E^,  because  the  line  integral  is  taken  over 
2rr  radians/ turn  of  the  probe  coil.  Therefore,  if  we  expand  E^  in  a 
Fourier  series,  {cos  n<{>},  then  only  the  n=0  term  will  contribute  a  non¬ 
zero  value  to  the  EMF  integral.  But  if  is  independent  of  <|>,  then  the 
E  vector  field  (which,  by  assumption,  consists  of  only  the  <t>-component) 
is  divergenceless.  This  is  equivalent  to  saying  that  the  Green's  function 
of  (6)  Is  the  electric  field  that  is  produced  by  a  circular  filament  of 
current,  of  radius  r' ,  located  at  the  plane  z=z' .  In  addition,  we  can 

also  take  the  "anomalous  current,"  J  ,  to  be  wholly  in  the  ^-direction, 

a 

and  divergenceless. 

Thus,  we  can  reduce  (6)  to  the  scalar  equation 


-  Ef)(r,z)  =  -jwyQ  ||  G(r,z;r' ,z')Ef(r' ,z') 


(aQ  -  of)(r',8')r'dr’dz'  , 


where  Eq,  E^,  and  G  are  the  ^-components  of  their  respective  fields.  The 
gist  of  the  preceding  argument  is  that,  by  using  coils,  as  shown  in 
Figure  2,  we  are  unable  to  determine  the  azimuthal  location  or  extent  of 
the  flaw.  In  carrying  out  the  integral  over  <J>  in  (6),  thereby  transforming 
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it  into  the  two-dimensional  integral,  (7),  we  assume  that  the  flaw  is 
confined  to  the  plane  4>=0.  Tuis  means  that  the  flaw  has  the  functional 
dependence 

(o  -  Of)(r,z,4>)  =  (0q  -  Cf)(r,z)  . 

In  oidei:  to  further  develop  the  model,  we  must  take  into  account 
that  the  region  of  interest  consists  of  three  parts:  the  interior  o£ 
the  tube,  the  tube  wall,  and  the  region  exterior  to  the  tube.  We  call 
these  regions  1,  2,  3,  respectively,  and  introduce  the  following  nota¬ 
tion  for  the  Green's  function: 

G^(r,2,  ;r' ,z')  =  Field  produced  at  (r,z)  in  region  i, 

due  to  filamentary  current  loop  at  (r'.z’) 
in  region  j 

=  G^Cr'.z’jr.z),  i,j  =  1,2,3  . 

The  last  equation,  a  reciprocity  relation,  follows  because  the  mag¬ 
netic  permeability  of  all  three  regions  is  the  same  (this,  and  other 
matters  relating  to  Green's  functions  in  stratified  media,  can  be  found 
in  [23]). 

We  introduce  the  following  notation 

E  7(r,z)  =  Electric  field  in  region  1  or  2,  with  flaw  present 

Eg(r,z)  =  Electric  field  with  flaw  absent,  due  to  exciting  coil. 

Then  (7)  produces  the  following  basic  integral  equation  for  computing 
E  in  the  flawed  region,  which  is  in  region  2: 

ff  °f 

E2(r,z)  +  jwp0a0  G22(r,z;r',z')E2(r',z')(^=-(r',z')  -  Dr'dr'dz' 
Flaw  0 

=  EQ(r,z)  (8) 
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In  addition,  we  have  the  integral  relation  for  computing  the  perturbed 
electric  field  at  the  probe  coil  (which  lies  within  region  1) : 

(Eq  -  E1)(r,z)  =  jiriUgOg  ||  G12(r,z;r,,z!)E2(r,,z')  • 

Flaw 

°f 

*  (~(r*  ,z')  -  Dr'dr'dz’  .  (9) 

When  this  equation  is  integrated  over  the  probe  coil  we  get  the  per¬ 
turbed  EMF.  If  we  assume  that  the  probe  coil  is  uniformly  and  densely 
wound  with  n^  turns  per  unit  area  (in  the  r-z  plane),  we  get  for  this 
EMF: 

r  f 

EMF  =  -2imc  jj  (Eg  -  E^rdrdz  .  (10 

Probe  Coil 

Finally,  the  electric  field.  Eg,  that  is  produced  by  the  exciting 
coil  is  given  by 

E0(r,z)  =  -ju)Pg2TT  |  ,z')Jg(r'  .z'Jr’dr'dz' 

Exciting 

Coil 

=  -jo)Pg2TineIg  ||  ,z;r'  jz'ir’dr'dz*  ,  (11 

Exciting 

Coil 

where  Jg  is  the  exciting  coil  current  density,  ne  is  the  density  of  tarns 
in  the  exciting  coil,  and  Iq  is  the  current  carried  by  the  exciting  coil. 

Equations  (8)-(ll)  constitute  the  model  system.  The  algorithm  for 
using  the  system  consists  of  first  computing  the  incident  field.  Eg,  at 
the  flaw,  by  (11) :  this  is  the  right-hand  side  of  (8) .  For  a  given  dis¬ 
tribution  of  flaw  conductivity,  a^(r,z),  (8)  can  be  solved  numerically. 

Its  solution,  the  electric  field,  E2>  in  the  flawed  region  is  the  source 
term  for  (9),  which  produces  the  perturbed  electric  field  at  the  probe 
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coil  in  region  1.  The  integral  of  the  perturbed  electric  field  produces 
the  perturbed  EMF,  (10),  which  is  then  compared  with  the  measured  EMF 
to  determine  if  the  assumed  flaw  conductivity,  a^(r,z),  is  "close"  to  the 
actual  (though  unknown)  flaw  conductivity.  The  problem  is  really  non¬ 
linear  because  (8)  involves  the  product  of  two  unknowns,  0^(r,z)  and 
E2(r,z).  Thus,  some  form  of  iteration  is  required,  in  which  one  starts 
with  an  assumed  distribution  for  a^(r,z),  and  then  hopes  to  converge  to 
a  final  acceptable  value. 

The  model  that  we  have  developed  is  quite  general.  Crucial  its 
application  as  a  reconstruction  technique  is  the  ability  to  compute  the 
Green's  function  of  the  "known"  region,  i.e.,  the  region  that  exists  in 
the  absence  of  any  flaws.  When  the  known  region  consists  of  a  cylindrical 
tube,  then  the  Green's  function  can  be  computed  in  a  straightforward  manner 
by  the  use  of  Fourier  transforms  and  algebra;  we  carry  out  these  computa¬ 
tions  in  Appendix  A. 

If,  however,  we  wish  to  reconstruct  flaws  that  exist  in  the  presence 
of  tube  supports,  as  in  Figure  3,  or  tube  flaring,  as  in  Figure  4,  or 
any  other  known  irregularities,  then  we  cannot  hope  to  compute  the 
Green's  function  in  a  purely  analytical  manner.  In  this  case  the  Green's 
function  satisfies  an  integral  equation,  which  must  be  solved  numerically. 
We  hardly  consider  this  to  be  a  very  serious  drawback,  however,  because 
so  much  of  our  modeling  effort  involves  the  numerical  solution  of  integral 
equations. 

The  integral  equation  that  is  satisfied  by  the  Green's  function  is 
identical  to  (8),  with  the  following  changes:  G^2  is  the  Green's  function 
of  the  cylindrical  tube,  as  if  there  were  no  irregularities  (the  sub¬ 
scripts  may  be  different,  depending  upon  the  location  and  type  of  known 
irregularity),  a^(r',z')  becomes  the  conductivity  of  the  known  irregular¬ 
ity,  and  the  integration  is  over  the  volume  occupied  by  the  irregularity, 
and  the  right-hand  side  is  replaced  by  G22*  Of  course,  we  are  still 
interested  in  only  the  <j>-component  of  the  Green's  function,  and  that  is 
why  we  can  use  the  scalar  integral  equation,  (8).  If,  in  (8),  there  is 
no  irregularity,  then  0^  =  a^,  and  the  Green's  function  is  identical  to 
that  for  the  circular  cylinder. 


t 
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An  additional  complication  exists  when  the  region  is  more  irregular 
than  a  cylinder,  and  that  is  the  need  to  approximate  the  region  in  order 
to  apply  the  method  of  moments  (or  any  other  numerical  method) .  In 
Figure  5  we  show  two  methods  of  approximating  a  curved  surface  with 
splines  of  zero-order  and  first-order.  Higher-order  approximating  func¬ 
tions  could  also  be  used,  but  the  two  shown  in  Figure  5  are  the  easiest 
to  work  with. 

(b)  Discretization  of  the  Model:  The  Method  of  Moments  [24] 

The  discretization  of  the  problem  via  the  method  of  moments  is  based 
on  the  use  of  a  mesh,  as  shown  in  Figure  2.  In  order  to  reduce  (8)  to  an 
algebraic  system,  we  expand  E^Cr.z)  and  (c^/Oq  -  1)  in  pulse  functions 
that  are  defined  with  respect  to  this  mesh: 


E-,(r,z)  ■  E  E  P.(r,z) 
j=l  3  3 


(12) (a) 


(~  (r,z)  -  1)  =  E  a.P.(r.z) 

°0  j-1  3  3 

where  N  is  equal  to  the  number  of  cells  in  the  mesh,  and  P  (r,z)  is  the 

£  J 

jth  pulse  function,  which  is  defined  by 


P^(r,z)  =  1  , 

=  0  , 


(r,z)  in  jth  ceil 
otherwise 


The  jth  expansion  coefficients,  E..,  a^,  are  the  constant  values  of 
the  fields  over  the  jth  cell- 

Because  E^  and  (a^/c^  ~  1)  have  identical  expansions  in  non-overlapping 
pulse  functions,  it  follows  that  their  product  does  also: 


E2(r,z)(of(r,z)/u0  -  1)  *  E  E^P^r.z)  . 
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Upon  substituting  (12)  and  (14)  into  (8),  we  get 


N  N 

c  c 

Z  E  P  (r,z)  +  jupna  E  E  a 

3=1  3  3  3=1  3 


ill 


G22(r,z;  r’,z')P  (r' .z'Jr’dr'dz' 


Flaw 


=  EQ(r, z) 


(15) 


Next,  we  take  moments  of  (15);  i.e.,  we  multiply  (15)  by  weighting 
functions,  Q^(r,z),  i  =  1,  .  .  .  ,  N^,  and  integrate  over  the  flaw.  If 
the  weighting  functions  are  delta  functions  that  are  located  at  the  center 
of  each  cell,  then  the  method  is  called  point -matching;  if  they  are  the 
same  functions  that  were  used  in  the  expansions,  (12),  then  the  method  is 
called  Galerkin's  method.  In  any  case,  the  result  is  the  svstem  of  N 

c 

equations: 


N 


c 


E  E 

J-l 


3 


Flaw 


(r,z)Qi(r,z)rdrdz 


N 

c 

+  E 

3 


(f 


aWVo  JJ 

Flaw 


rdrdz 


Flaw 


G22(r,z;  r’ ,z')Pj (r* ,z’)Q^(r,z)r,dr,dz' 


|  EQ(r,z)Qi(r,z)rdrdz  , 
Flaw 


(16) 


The  vector-matrix  version  of  (16)  is: 
(A  +  3w}i0a0G)E  =  F  , 

where 

Aij  *  {}  z) rdrdz 

Flaw 


(17) 


(18) 
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|  Qi(r,z)rdrd2  Jj  G22(r,z;r'z')Pj  (r'.z^r'dr'dz' 
Flaw  Flaw 

Fi  *  j{  E0(r,z)Qi(r,z)rdrdz  ,  i,  j  »  1,  .  .  .  ,  H£ 


Our  analysis  has  been  based  on  the  method  of  point-matching,  with 
Qi(r,z)  *  6(r  -  r.)S(z  -  z^/r,  where  (r^ , z^)  are  the  coordinates  of  the 
midpoint  of  the  ith  cell.  Point -matching  generally  takes  better  advantage 
of  the  singularity  in  the  Green's  function  to  produce  a  better  conditioned 
matrix  (i.e.,  one  more  diagonally  dominant)  for  inversion.  The  disadvan¬ 
tage  of  point-matching  is  that  the  infinite  integrals  that  define  the 
various  matrix  elements  in  G_  and  F^  (see  Appendix  B)  do  not  converge 
as  rapidly  as  with  Galerkin's  method,  so  one  has  to  take  more  care  in 
their  numerical  evaluation.  This  has  not  turned  out  to  be  a  problem, 
however. 

“I* 

Thus,  upon  letting  (z_, ,  z.)  be  the  lower  and  upper  z-limits,  re- 
s  J  J 

spectively,  and  (r^ ,r^)  the  lower  and  upper  r-limits,  respectively,  of 
the  jth  cell,  we  get 


+  + 
r*  Z\ 

j  Jj  3  <S(r-r^)6(z-z^) 

Vi 

+  + 

rj  Z1 


rdrdz  *  6, 


.  (19) 


eu*‘,J  I  G22<Wc'-2’>l’d'’d2' 


r.  z. 
3  3 


Fi  *  E0(ri»V  =  If  G21(ri,z^;r' .z^r'dr'dz' 


Exciting 

Coil 


where  5^  *  1,  if  i  -  j,  and  -  0,  if  i  #  j.  Similarly,  when  (14)  is 
substituted  into  (9),  we  get: 


I 


» 


i 


i 


* 


% 


I 


+  + 

Nc  rj  zi 

(Eg  •  E^) (r,z)  »  j^PgOg  2  J  j  G12(r,z;r* .z'Jr’dr'dz*  ,  (20) 

r. 

J  j 

Now  we  define  the  kth  probe  coil  as  one  having  an  inner  radius,  p^, 

outer  radius,  p„,  left-hand  z-coordinate,  £,  right-hand  z-coordinate, 

+  *  K 

Z^,  and  midpoint  z-coordinate  of  5^.  Then  the  EMF  induced  into  the  kth 
probe  coil  is  given  by  substituting  (20)  into  (10) : 


Nc  v2  ^k  ‘j  ‘j 

-j(oy0a02imc  2  |  rdrj  dz  J  j  G12(r,z:r’ .zOr’dr'dz’ 


+  + 

Vi 


EMF(k) 


P1  ^k 


rJ  *3 


k  *  1,  ■  *  •  ,  K  • 


(21)- 


This,  too,  can  be  put  into  vector-matrix  form: 

EMF  «  T (OE)  ,  (22 % 

where 

„+  +  + 
p2  ?k  rj  zj 

Ty  *-ju)ygag2imc|rdr|  dz  |  j  G12(r,z;r’ .z^r’dr'dz*  (23) 

P1  ?k  rj  zj 

is  the  transfer  function  from  the  jth  cell  to  the  kth  coil.  Note  that 

the  number  of  coils,  K,  is  not  necessarily  equal  to  the  number  of  cells, 

N  .  Indeed,  in  applying  the  method  of  least-squares  to  this  model,  we 
c  - 

take  K  *  50,  which  produces  100  real  and  imaginary  components  of  EMF* 
and  N£  *  60.  This  yields  an  overdetermined  system,  which  is  typical  of 
least-squares  problems. 

In  Appendix  B  we  derive  egressions  for  the  various  vector-matrix  ;  . 
elements  that  have  just  been  defined.  ~  ‘  J  -X 
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III.  LEAST-SQUARES  ALGORITHMS:  UNPACK  AND  MINPACK 


Now  that  we  have  a  model  system  in  hand,  we  shall  say  more  about  the 
mathematical  and  numerical  solutions  of  linear  and  nonlinear  least-squares 
problems.  We  rewrite  (1)  as  the  vector  equation 

e  =  f  (a)  (24) 

We  seek  a  least-squares  solution,  as  defined  below.  Let 

F(c)  =  f(0)  -  i  .  (25) 

Upon  introducing  the  usual  squared-norm  notation,  we  have 

I  |F(o)|  1 2  =  ||£(o)  -  e||2  -  [  Z  [f^,  ...  ,  cM)  -  ei]2]1/2  .  (26) 

Then  the  definition  of  a  least-squares  solution  of  (24)  is:  given 
e  =  (e^,  ....  eN),  find  o  *  (01,  .  .  .  ,  0M)  that  minimizes  ||f||2; 
i.e.,  solve 

min| |f (5)  —  e | | ^  •  (27) 

0 

We  first  consider  the  case  where  f  is  linear.  Then  write  (24), 

(25),  and  (27),  respectively,  as 

e  =  A  0  (28) 

r  *  e  -  A  5  (29) 


min|  |e  -  A  <j|  L 
0 


If  5  is  a  solution  of  (30),  then  it  is  known  that  [25] 


aT  — *  *T  —  *— * 

Aa  r  =  A  (e  r  A0  )  *  0  , 
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where  the  superscript,  T,  denotes  the  transpose  of  a  matrix.  Thus 


=T=-*  =T- 

A  AO  =  A  e 


(32) 


or 


-* 

o 


=+  - 
=  A  e 


(33) 


s*4»  =Tsa  —  S 

where  A  =  (A  A)  A  is  the  pseudoinverse  of  A.  While  (33)  characterizes 
the  mathematical  solution  of  (30),  we  don't  actually  numerically  compute 
A  ;  for  numerical  solutions  other  methods  are  used. 

A  standard  method  of  solving  (30)  is  to  use  the  QR-factorization  of 
the  matrix  A.  Given  that  A  is  N  x  M  (where  N  M),  there  exists  an 
orthogonal  matrix  Q  (of  order  N  x  N)  such  that 


(34) 


where  R  is  upper  triangular.  If  we  write  Q  =  ,  where  has  M 

columns ,  then 


A  =  QXR 


(35) 


Thus  if  rank  (A)  =  M,  the  columns  of  Q  form  an  orthonormal  basis 

s s  s  as  ss 

for  the  column  space  of  A.  Now  if  A  *  [A,,A„],  where  A  has  k  columns, 

1  L.  1 


and  if  R 


*11 
0  R, 
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22 


where  R^  is  kxk,  then 


(36) 


Hence  Q  and  R^  give  a  QR-factorization  of  A^.  This  truncated 
decomposition  is  important  for  matrices  whose  rank  is  less  than  full, 

a 

i.e.,  for  which  rank  (A)  <  min(M,N). 
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This  orthogonal  triangularization  is  generated  by  Householder  trans¬ 
formations  [21,  25].  Once  we  have  this  triangularization,  then  the  solu¬ 
tion  to  (30)  follows: 


=T—  =T-  =T=- 

Q  r  =  Q  e  -  Q  AO 


sT-  1 
Qr. 


R  „  .  1  =T- 

—  a,  where  -  =  Q  e 

0  e2 


=T-  ei  “ 

Qr.  ^ 


Since  Q  is  orthogonal,  |  j Q  ir |  1 2  =  llrj^.  Thus  |  1  r |  1 2  is  minimized 


-  Ro  , 


and  when  this  is  true 


I  e2 1  ^  2 


It  is  important  to  realize  that  when  k  =  rank(A)  =  min(M,N),  then 
the  two  subroutines  DQRDC  and  DQRSL  from  LINPACK  [21]  provide  the  method 
for  solving  (30) .  DQRDC  produces  the  QR  decomposition  of  A  with  column 
pivoting.  This  resulting  matrix  is  passed  to  DQRSL  for  the  solving 
stage.  DQRSL  uses  the  k  =  rank(A)  columns  to  produce  the  desired  least- 
squares  fit. 

If,  however,  A  is  rank  deficient  (or  near  rank  deficient),  then  we 
need  a  truncated  least-squares  fit.  This  can  be  achieved  by  using  the 
subroutine  DQRST  [21,  page  9.11].  This  subroutine  allows  for  a  user 
supplied  tolerance  and  calls  DQRDC  and  DQRSL.  Based  on  this  tolerance, 
some  of  the  columns  of  the  output  of  DQRDC  are  zeroed.  DQRSL  then  pro¬ 
duces  the  truncated  least-squares  fit. 


1.,  hpi  nj  ;i  ii.'jnii  111  ']i<  iljlli  i]j  'j[!l!^ill  ‘I'ii  j  i|i  H!ii;ini,lilll,ilijlii!iii!i)  1,1}  ill  iifl  ilitlljltl]  I  liilljiJ  'liiijlL  Ifa  1  'Jiil  Ife 


We  apply  the  linearized  theory  that  has  just  been  described  to  (22), 
with  the  electric  field  column  vector,  E,  replaced  by  its  zeroth  order 
approximation,  E^,  so  that  the  resulting  equation  becomes 

EMF  =  (T,  .E  ,)o  .  (42) 

k  k.]  oj  j 

When  we  generated  a  singular- value  decomposition  of  T,  we  found  that 

the  ratio  of  the  largest  to  smallest  singular  values  was  on  the  order  of 
12 

10  .  This  ratio  is  a  measure  of  the  condition  of  the  matrix.  Thus, 

while  we  were  working  with  a  highly  ill-conditioned  matrix,  our  model 
produced  very  accurate  results  for  inversion,  as  we  will  show  in  the  next 
section. 

We  have  pointed  out  before  that  cur  model  equations,  (8)-(ll),  or 
their  discretized  versions,  (17)-(23),  are  in  reality  nonlinear.  Thus, 
we  must  use  a  nonlinear  least-squares  scheme.  The  subroutines  LMDER  and 
LMSTR  in  MINPACK  [22]  fill  this  bill  nicely.  Both  are  based  on  the 
Levenberg-Marquardt  method,  which  we  briefly  describe  below  [26,  27]. 
Recall  our  problem:  (25)-(27).  If  we  linearize,  we  see  that 


1  |F(cr  +  £)||2  =  ||F(5)  +  F'  (5)5||  2  -  *(P),  (43) 

where  F*  (o)  is  the  Jacobian  matrix.  A  standard  way  of  minimizing  i[i(p) 
is  by  using  a  Gauss-Newtown  method: 


given: 
solve: 
then  let: 


0  (k),F(0^k)),J(a(k))  =  F'(0(k)) 

J(5^)p^  =  -F(o^)  in  a  least-squares  sense 


-(k+1)  .  -(k)  +  -(k) 


This  works  well  if  J(a)  is  full  rank.  In  the  rank  deficient,  or 
nearly  rank  deficient  case,  however,  modifications  are  required.  It 
should  be  noted  that  in  practice  rank  deficiency  arises  often. 

Since  the  linearization  of  (43)  is  not  valid  globally,  we  consider 
a  constrained  linear  least-squares  problem: 
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:  |  |  Dp  |  1 2  _<  A}  »  min{  |  |  F  +  Jp  |  j  2 :  |  |  Dp  |  1 2  £  A)  , 

P  P 

=* 

where  D  is  diagonal.  The  Levenberg-Marquardt  method  is  based  on  the 

— Jjl  — jjjj  (— 

fact  [26]  that  if  p  i3  a  solution  to  (44),  then  p  =  p(X),  for  some 
X  ^  0,  where: 

=Ts  =T=  -  =T- 

(J  J  +  XD  D)p  =  -J?  . 

The  way  to  solve  for  p  is  to  recognize  that  (45)  are  the  normal 
equations  for  the  linear  least-squares  problem 


J  -  _  _  F 

*1/2d  P'  5 


The  implementation  of  these  facts  is  the  basis  of  the  nonlinear 
least-squares  subroutines  in  MINPACK. 

In  order  to  apply  the  nonlinear  least-squares  algorithm  to  our 
discretized  model,  we  must  send  to  LMDER  and  LMSTR  the  nonlinear  func¬ 
tion  to  be  minimized,  as  well  as  the  Jacobian.  The  Jacobian  is  obtained 
from  (22)  as 


J«"Vj  • 


which  is  the  matrix  that  multiplies  the  column  vector  a.  Therefore, 
the  steps  to  be  followed  in  solving  the  nonlinear  least-squares  problem 
are: 

•  enter  initial  guess  for  0  and  solve  (17)  for  E 

•  compute  Jacobian  from  (47) 


•  compute  EMF^el  ^rom  (^2) 


form  F(a)  »  EMF  .  .  -  EMF  , 

model  measured 


*  call  LMDER  or  LMSTR 
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IV.  EXAMPLES  OF  RECONSTRUCTION  OF  SIMULATED  FLAWS 

The  theory  of  inversion  involves  two  components,  a  theoretical  model 
that  is  based  on  a  rigorous  application  cf  electromagnetic  theory,  and 
numerical  algorithms  that  effectively  implement  least-squares  theory. 

Each  of  these  has  been  dealt  with,  and  now  we  illustrate  how  the  method 
works  for  the  reconstruction  of  computer  simulated  flaws. 

All  numerical  experiments  were  run  in  double  precision  on  the  PRIME 
550-11  and  IBM  360  machines.  The  double  precision  data  word  on  the  PRIME 
occupies  64  bits,  of  which  47  are  the  mantissa  and  16  the  exponent.  The 
effective  precision  is  about  14  digits.  The  IBM  double  precision  word 
has  a  56  bit  mantissa,  which  allows  an  effective  precision  of  about  17 
digits.  Precisions  such  as  these  are  required  for  meaningful  computa¬ 
tions,  because  the  condition  number  of  the  Jacobian  matrix  is  phenomenal — 
12 

on  the  order  of  10  .  Even  with  this  large  condition  number,  the  compu¬ 

tations  produced  excellent  results;  in  the  worst  case  the  reconstructions 
were  exact  to  at  least  three  places  on  the  PRIME,  and  five  places  on  the 
IBM.  This  verifies  that  the  algorithms  in  the  UNPACK  and  hINPACK  pack¬ 
ages  tend  to  work  better  in  higher  precision. 

The  physical  system  that  was  modeled  is  a  variation  of  the  multi¬ 
coil  system  that  has  been  described  earlier,  and  is  shown  in  Figure  6. 

It  consists  of  a  fixed  exciting  coil  and  a  single  probe  coil  that  can 
be  moved  axially.  This  system  is  typical  of  a  common  flaw  detection 
scheme-  The  mesh  on  which  the  discretization  is  defined  is  also  shown. 

It  consists  of  six  rows  of  ten  cells,  and  spans  the  entire  tube  wall- 
thickness.  The  starting  position  of  the  probe  coil  is  at  the  left  edge 
of  the  mesh,  and  the  final  position  is  at  the  right  edge.  The  probe 
coil  is  stepped  through  fifty  equal  intervals  between  these  limits,  . 
thereby  generating  a  total  of  100  real  and  imaginary  EMF  values  that  are 
used  in  the  least-squares  inversion. 

The  physical  parameters  of  the  model  are  typical  of  real  systems. 

The  inner  radius  of  the  tube  is  0.510",  and  the  outer  radius,  0.375". 

The  length  of  the  mesh  is  0.50"  in  the  z-direction,  thereby  giving  a 
cell  resolution  of  0.05"  by  0.011".  The  probe  coil's  inner  radius  is 
0.05",  outer  radius,  0.100",  and  its  length  is  0.50",  the  same  as  the 
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mesh.  The  exciting  coil  is  centered  on  the  mesh  in  the  2-direction 
(neither  of  these  last  two  items  is  a  requirement  of  the  inversion 
method).  The  density  of  turns  of  the  exciting  coil  is  2x10^  tums/m, 
which  is  comparable  to  that  of  20  gauge  copper  wire.  The  probe  coil 

has  an  inner  radius  of  0.100",  outer  radius  of  0.26",  and  a  length  of 

0.250".  Its  turns  density  is  2x10^  tums/m,  which  is  comparable  to  30 
gauge  copper  wire.  The  tube  conductivity  is  3.5x10^,  which  is  equal 
to  the  conductivity  of  alumninum,  and  the  freuqency  of  operation  is 
1kHz. 

In  Figure  7  we  show  a  simulated  flaw  (the  "original")  at  the  top  and 
its  reconstructed  version  at  the  bottom.  The  real  (R)  and  imaginary  (I) 
parts  of  the  perturbed  EMF,  as  measur  >d  by  the  probe  coil  when  it  is 

moved  acorss  the  mesh,  are  shown  in  the  middle  of  the  figure.  This  EMF 

curve  is  actually  an  interpolation  based  on  the  fifty  probe  coil  positions. 
In  this  figure,  and  the  next  two,  we  simulate  the  flaw  by  letting  =>  0 
at  the  flaw  location,  and  0^  *  off  of  the  flaw.  Thus,  according  to 


(12) (b),  a. 


-1  if  the  jth  cell  lies  on  the  flaw,  and  =  0,  otherwise. 


Note,  in  Figure  7,  that  because  the  original  flaw  is  placed  sym¬ 
metrically  in  the  mesh,  the  EMF  is  symmetrical  about  the  center  of  the 
mesh,  also.  The  reconstruction  is  clearly  perfect  (to  at  least  three 
significant  digits) ,  indicating  that  the  least-squares  inversion  algo¬ 
rithms  work  quite  well  in  this  model.  We  must  be  cireful  to  note,  how¬ 
ever,  that  in  this  report  we  have  considered  only  original  flaws  that 
are  defined  on  the  same  mesh  as  that  used  for  reconstruction;  i.e.,  each 
part  of  the  flaw  is  constant  over  a  full  cell  of  the  reconstruction  mesh. 
We  intend  to  consider  the  more  general  case,  in  which  the  flaw  may  be 
defined  on  a  different  mesh  than  that  used  for  reconstruction  (say,  one 
with  smaller  cells,  or  cells  that  are  displaced  from  the  cells  of  the 
reconstruction  mesh) .  This  will  test  the  ability  of  the  model  to  re¬ 
solve,  as  well  as  invert,  data. 

To  satisfy  ourselves  that  the  excellent  results  that  were  obtained 
in  Figure  7  were  not  due  to  symmetry,  we  considered  the  asymmetrical 
flaws  of  Figures  8  and  9.  Again,  the  reconstruction  was  perfect  to  at 
least  three  signficicant  digits.  It  should  be  noted  from  these  three 
examples  that  the  more  concentrated  the  flaw,  the  greater  is  the  peak 
of  the  EMF  curve. 
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A  crucial  test  of  inversion  in  highly  ill-conditioned  systems  has 
to  do  with  corrupted  data.  The  question  is,  does  the  reconstruction 
"follow"  the  corrupted  data,  or  does  the  result  lose  all  significant 
figures?  In  order  co  test  our  model’s  response  to  corrupted  data,  we 
performed  the  following  numerical  experiment.  We  assigned  to  each  cell 
in  the  mesh  a  number  between  0  and  1,  chosen  at  random  by  using  the 
FORTRAN  random  number  generator.  Then  the  model  EMF  that  is  produced  by 
this  "flaw"  is  computed  by  following  the  first  three  steps  that  are 
listed  below  (47) .  This  "true"  data  is  then  corrupted  by  adding  to  it 
the  same  data  multiplied  by  either  0.01,  0.10,  or  0.20,  and  then  using 
this  as  the  "measured"  EMF  (see  the  discussion  following  (47)).  Figure 
10  shows  the  results  of  this  experiment.  There  we  show  the  original 
flaw,  consisting  of  the  sixty  randomly  chosen  cell  conductivities, 
followed  by  the  reconstructed  flaw,  simulated  by  the  sixty  values  of 
computed  cell  conductivities,  for  the  case  of  1%,  10%,  and  20%  corrupted 
data. 

Again,  the  results  are  excellent.  We  don't,  of  course,  expect  to 
reconstruct  the  original  flaw  by  using  corrupted  EMF  data.  We  are 
happy,  though,  to  see  that  the  reconstructed  flaw  "tracks"  the  original 
flaw,  in  the  sense  that  it  departs  by  almost  exactly  1%,  10%,  or  20% 
from  the  original.  Such  stability  in  the  face  of  a  very  ill-conditioned 
system  attests  to  the  excellence  of  the  UNPACK  and  MINPACK  algorithms. 

The  same  results  that  are  shown  in  Figures  7-10,  are  obtained  with 
either  the  linear  or  nonlinear  algorithms  that  are  described  in  Section 
III.  The  reason  for  this  is  that  in  (17)  the  term  involving  the  matrix 
G  is  much  smaller  than  the  first  term,  A.  Thus,  the  solution  of  the 
equation  is  £  :  Eq,  and  when  this  is  substituted  into  (22),  or  (47),  we 
see  that  the  Jacobian  matrix  is  constant,  so  that  the  nonlinear  algo¬ 
rithms  may  be  replaced  by  the  simpler  linear  ones. 

V.  COMMENTS  AND  CONCLUSIONS 

In  this  report  we  have  developed  a  model  for  eddy-current  inversion 
that  is  based  on  the  application  of  rigorous  electromagnetic  theory  and 
numerical  algorithms  for  least-squares.  So  far  we  have  attempted  only 
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to  verify  the  inversion  method  for  flaws  that  were  defined  on  the  same 
grid  as  that  used  for  reconstruction.  The  results  have  been  excellent, 
and  suggest  that  the  method  can  be  used  as  the  basis  for  the  development 
of  an  engineering  prototype  system.  Before  such  a  system  can  be  ef¬ 
fected,  however,  we  believe  that  the  following  additional  studies  should 
be  carried  out: 

*  study  resolution  of  flaws  that  are  defined  on  a 
different  grid  than  that  used  for  reconstruction 

*  carry  out  examples  of  reconstruction  in  the  presence 
of  known  irregularities 

*  determine  computer  hardware  requirements 

*  determine  computer  software  requirements 

*  determine  relative  advantages  of  multicoil,  multi¬ 
frequency,  and  transient  (time-domain)  systems 


optimize  exciting  coil  and  probe  configurations  and 
design. 
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Figure  1.  Illustrating  three  systems  problems:  (a)  the  "direct" 
problem,  (b)  the  "signal-detection"  problem,  and 
(c)  the  "inverse"  system-identification  problem. 


other  known  irregularity. 


Figure  5.  Illustrating  the  approximation  of 
a  curved  surface  by  splines  of 
(a)  zero  order  and  (b)  first  order 


EXCITING  COIL  IFIXEDI 


exciting  coil  and  a  single  orobe  coil  that 
can  be  moved  axially  across  the  mesh. 


RECONSTRUCTED  FLAW 

Figure  7.  Illustrating  a  symmetrically  placed  flaw  (top),  the  real  (R) 
and  imaginary  (I)  parts  of  the  EMF  induced  into  the  probe 
coil  (center),  and  the  reconstructed  flaw  (bottom).  The  flaw 
consists  of  the  darkened  cells.  The  reconstruction  is  exact 
to  the  least  three  significant  digits. 
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Figure  9.  Illustrating  the.  reconstruction  of  another  asymmetrically 
placed  flaw. 
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Figure  10.  Illustrating  the  reconstruction  of  random  flaws  with 
perturbed  EMF  data:  (1)  1%  perturbation,  (b)  10% 
perturbation,  (c)  20%  perturbation. 


APPENDIX  A.  CALCULATION  OF  THE  GREEN'S  FUNCTION 


The  interest  of  the  cylinder  is  labeled  region  1,  the  tube  wall, 
region  2,  and  the  exterior  to  the  cylinder,  region  3.  In  each  region 
the  magnetic  permeability  is  Pq  and  the  dielectric  constant,  £q.  The 
electrical  conductivity  of  regions  1  and  3  is  zero,  whereas  that  for 
region  2  is  0.  We  use  the  notion 


G„(r,z;r'  ,z')  =  field  produced  at  (r,z)  in  region  i. 

due  to  a  filamentary  current  loop  at 
(r',z')  in  region  j,  where  i,j  =  1,2,3. 

Using  this  notation,  the  Green's  function  satisfies 
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If  we  write  G , .  =  G.  .a,,  where  the  scalar  G.^  is  independent  of  d>, 
ij.  ij  <P  ij 

because  the  vector  G^.  is  divergenceless,  then  the  vector  differential 
operator  on  the  left-hand  side  of  these  equations  is  replaced  by  the 
scalar  differential  operator 


3g. 


3r 


ij 


(A-4) 


where  k  stands  for  either  or  k2- 

The  only  parts  of  the  Green's  function  that  are  required  in  the 
model  equations,  (8)-(ll),  are  G^*  G^>  30(1  G^;  but  G21(r' ,z*  ;r,2>  = 
G^2(r, z;r' , zr)  because  M  »  in  each  of  the  three  regions  (see  [23]  of 
the  main  text).  Hence,  all  we  need  be  concerned  with  is  (A-2) ,  which 
we  rewrite  using  (A-4) : 
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We  use  Fourier  transforms  to  solve  (A-5).  Thus, 
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so  that  (A-5)  becomes 
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where  we  have  used  the  completeness  relation 


-^f  e-Jh<z-'’)dh  , 


(A-8) 


and  ^^*1,  ifi«2,  and  =  0,  otherwise. 

Since  (A-7)  holds  for  all,  2,  z',  we  can  equate  integrands  to  get 
the  differential  equation  that  is  satisfied  by  the  transformed  variable: 
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I  et  r^  be  the  inner  radius  of  the  cylinder  and  r^,  the  outer.  Then 
the  boundary  conditions  that  are  to  be  satisfied  by  the  G are: 


n  =>  1,2 


(A- 10)  (a) 
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Equations  (A-10) (a) , (b)  imply  the  continuity  of  tangential  electric 
and  magnetic  field,  respectively,  across  a  surface  that  does  not  carry 
magnetic  or  electric  current  singularities.  Equation  (A-10)(d),  on  the 
other  hand,  implies  that  the  tangential  component  of  magnetic  field 
intensity  suffers  a  discontinuity  of  amount  -l/4iTr'  in  crossing  the 
surface  r  =  r',  on  which  the  filamentary  current  source  resides.  Alter¬ 
natively,  (A-10) (b) , (d)  can  be  obtained  by  integrating  (A-9)  an  infini¬ 
tesimal  distance  across  the  surfaces  r  «  r^,  r  *>  r^  or  r  *  r',  and  then 
invoking  the  continuity  of  electric  field,  as  expressed  in  (A-10) (a) , (c) , 
across  these  same  surfaces. 
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Because  the  source  point,  r  *  r',  is  in  region  2,  we  have  the  follow¬ 
ing  solutions  of  (A-9) : 


Gi2  -  AJ^Ctyr),  r  <  rL  ,  aQ  =  (k^  -  h2)1/2 


(A-ll) (a) 


G22  =  BH^C^r)  +  Ch£2)  (c^r) ,  ^  <  r  <  r’ ,  =■  (k2  -  h2)1/2  (b) 


G22  =  DH^L)  (c^r)  +  EH<2)(ctir),  r'  <  r  <  r2 


(c) 


532  '  »i2)<V)’ 


(d) 


where  A-F  are  arbitrary  constants  that  are  chosen  to  satisfy  the  boundary 
conditions  (A-10) . 

We  use  the  Bessel  function  in  region  1,  because  that  region  in¬ 
cludes  the  z- axis,  r  *  0,  and  J  is  regular  there.  The  Hankel  function, 

(2)  1  (2) 

H,  ,  is  used  in  region  3 ,  because  that  region  extends  to  r  «  00 ,  and  H 

1  (2)  x 

is  regular  there  (that  is  represents  an  outgoing  cylindrical  wave 

at  infinity). 

The  six  constants,  A-Fr  are  determined  by  applying  (A-10)  to  (A-ll), 
with  the  result: 

AJ1(aQr1)  =  BHj1)(a1r1)  +  CH<2)  (a^)  (A-12)(a) 


WVi>  ■  rf'Vl1  +“iCh'2)<Vi> 

SH^G^r')  +  CH^^^r’)  «  BH^^r’)  +  EH^^^r’) 
a1BHj1)(a1r*)  +  «1ChJ2)  (c^r*)  -  c^DhJ15  (a^’)  +  c^EhJ2*  (o^’)  + 


(b) 


(c) 


1 

2 

4tt  r’ 


(d) 
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DH^1)(a1r2)  +  EH<2)  (a^)  =  FH^2)  (aQr2) 


(e) 


“iDBi1)(V2)  +  “rfSV  ■  “orHo2)<V2)  ■ 


(£) 


In  arriving  at  the  final  form  of  these  equations,  we  used  the  Bessel 
function  identity  dB^/dz  +  B^/z  =  B^,  where  B  stands  for  J,  . 

These  equations  can  be  solved  in  a  straightforward  manner,  and  th 
constants  then  substituted  back  into  (A-ll) .  The  results  are: 


V.H^V.r’)  +  V9E(2)<a r') 

6.  „(r,h;r')  = ^ ± - J,  (a_r) 

U  (2*)  r  (V  V  -V  V  ) 


(A-13) (a) 


G22(r,h;r‘) 


[V^15  (ajr*  )+V,H<2)  (c^r ' )  ]  [V3H^1)  (air)+v4H^2)  (o^r)  ] 
-jlW^-V^) 


r  1  r  <  r’ 


(b) 


G22(r,h;r') 


[V3H^1)  (OjT ■ ' )+V4H^2)  (axr ' * ) ]  [VjH*3^  (c^rHv^p (c^r) ] 


-j/6TT(V1V4-V2V3) 


r'  <  r  <  r 
r  -  2 


(c) 


V  H(1)(a  r’)+V  H(2)(a  t’) 

S..2(r,hir')  --3-1-  ,-1 - 4-1 - 4—  SL^Vr)  , 

VVrtV 


(d) 


where 


Vi  "  alHl~  ^a0r2^H0  ^alr2^  “  a0H0  ^a0r2^Hl  *  ^a\r2^ 


(A-14) (a) 
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G22(r,2it',z') 


16ir 


Ev3H^1)  (ct^r*  )+V4h|2)  (OjV  )  ]  [V^15  (OjTH^hJ25  (c^r)  ] 


(V1V4~V2V3) 


3-jh(z-z’) 


dh 


r'  <  r  £ 


(c) 


G32(r,z;r'zf)  = 


-3-1  --v1- - ^ - i - Hj2)(aftr)e  jh(z  2  5dh 

<2*>  VVlVV2V3> 


(d) 


APPENDIX  B.  CALCULATION  OF  MATRIX  ELEMEN1S 

In  this  appendix,  we  derive  expressions  for  the  forcing-function 
vector,  (19) (c),  and  the  matrices,  (19) (b),  and  (23). 

In  order  to  compute  (19) (c) ,  we  start  with  the  reciprocity  relation 

C2l<r,I!r',z,>  •  G^r'.z'ir.z)  ,  CB-1) 

where  the  last  equality  follows  from  (A-15)(a).  We  assume  then  that  the 

(a\  (a)  (a)  (e\ 

exciting  coil  occupies  the  region  '  £  r  £  P2  ,  £  z  £  52  * 

Then,  when  we  substitute  (B-l)  into  (19) (c),  and  interchange  the  order  of 
integration,  we  get 


( 


-J^oVe 

S  ■— —  — 

TTr, 


§ 


§ 


§ 


€ 


* 


t  00 


.(e) 


e"jh(VC  )sin{hl,e/2)l(p^e)a0,p^e)a0) 


ha^ (v1v4-v2v3) 


•  [VjH^Ca^)  +  V2H^2)(a1ri)]dh  , 


(B-2) 


(e) 

where  Lg  is  the  length  and  £  the  midpoint  of  the  exciting  coil,  and 


■ 

Uz2,zx)  *  J^KdC  . 


(B-3) 


Before  computing  (19) (b) ,  we  introduce  the  vector  notation 


H,  -  V,H,(1)  +  VnH,(2) 


”12  “1  '  ’1“1  •  ’2“1  »  V3A  ’  H1  “  V3H1  *  T4“l 


(1>  +  V,H?> 


In  addition,  we  let  x>  be  the  larger  of  (r^.r') ,  where  r^  is  the 
radial  midpoint  of  the  ith  cell,  and  r<  the  smaller.  Using  this  notation 
allows  us  to  write  (19) (b)  concisely  as 


+  + 
r.  z. 

,1,3 


G  **  0 

«  j 


G22(ri’zi;r' .*,)*'dr,d** 


rj  Zj 


i 


I 


43 


t 


,Zi 


°  4-f  J  f 

j  16TT  J  '  J 

_<»  — 
•9 


r  '  iVVV«f34Vl^ 
e  jh(zi"z  }dz'  I  (V1V4'V2V3} 


®Je"jh{VZ1)-e'jh(VZi) 


(V12*^1  ^l*^  ^  ^34*^1  (ctir<)  ^  , 

(viV¥3) 


i  I  ^ 

f  F  ^ 


i  I  » 


In  order  to  evaluate  the  inner  integral  we  must  consider  three 

mm 

cases:  (i)  r^  <  r^,  (ii)  r^  -  r.,  and  (iii)  ^  >  r^  .  In  (i),  clearly 
r<  *  r^,  v>  »  r’,  and  in  (iii)  r<  *  r* r>  *  r^.  Hence,  the 
inner  integral  above  can  be  immediately  computed  for  these  two  cases: 


(V12-11l<V'»<V34-8l(Vl»  , 

(VlVV2  V  r 


(v12-  fl(vt.v;»(T34-H1(v<)) 

+ 

fr3  tvVVi^VV'” 

J  - <WVV -  r  ' 


(Vl2-Hi(airi))(V34-  H^r,,^)) 

ai(VlV4~V2V3) 


4  •  » 
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. . . . * . . . . . . 


The  vector  H  is  defined  by 


““’(Vi*  ■  f  Hf 


(B-5) (a) 


(z2,21) 


1 


(OM5 


For  case  (ii)  we  have  for  the  inner  integral 

+ 

, r±  (vvy^tys^v'))  t 
J.  (VrtV  1  (viVviV 


^V12*^l^airi^  ^34*  ^^airi,Ctlri^  +  ^12*  ^^airi*airi^  ^34*®l^“lri^ 


Hence,  when  these  results  are 


««  *  ^  S  f  Miai 


ij  J  Sit 


j  fl  |  "g-jhCz^z^)  sin(hA/2 


w 

j  °j  |  e-jh(zi-2^)  sin(hA/2) 


I 


III 


where  (i) ,  (ii) ,  (iii)  are  the  corresponding  functions  presented  above 
for  the  results  of  the  radial  integration,  and  A  is  the  width  of  a  cell 
in  the  z-direction.  We  point  out  that  the  mesh  with  which  we  are  working 
is  regular,  in  the  sense  that  its  ce!ls  are  of  constant  width  in  the 


r-  and  z-directions. 


The  expression  for  is  also  easily  derived  by  substituting 


(A-15)(a)  into  (23)  and  interchanging  the  orders  of  integration,  with 


the  result 


f  \2‘  ICP2cto»°ian) 

1  „  a!ct?(V  V  -V  V  ) 

0  0  ,  1  4  2  j 


kj  irr. 


cos  [hC^-z^)] 


cos[h(?^-z*)]  -  cos [h(t^-z”) ]  +  cos[h(?£-z*)] 
h2  • 


dh  (B-7) 
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